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Abstract 

The theory of a novel bond-order potential, which is based on the block Lanczos 
algorithm, is presented within an orthogonal tight-binding representation. The block 
scheme handles automatically the very different character of a and vr bonds by intro- 
ducing block elements, which produces rapid convergence of the energies and forces 
within insulators, semiconductors, metals, and molecules. The method gives the first 
convergent results for vacancies in semiconductors using a moments-based method 
with a low number of moments. Our use of the Lanczos basis simplifies the calcu- 
lations of the band energy and forces, which allows the application of the method 
to the molecular dynamics simulations of large systems. As an illustration of this 
convergent 0{N) method we apply the block bond-order potential to the large scale 
simulation of the deformation of a carbon nanotube. 
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1. INTRODUCTION 



To understand meso- and macro-scale phenomena from the atomistic level is an important 
subject in computer-aided materials modelling. This challenging study is not only intended 
as a realistic search for useful materials, but also for finding novel cooperative phenomena 
involving many atoms within large systems. Computer simulations of materials have in- 
evitably promoted the development of efficient algorithms for dealing with long time-scale 
phenomena. These methods have developed via two different approaches. The first, based 
on the continuum mechanics, is a hybrid approach that combines the continuum mechanics 
with atomistic simulations.0'i The second is more directly applying the molecular dynamics 
(MD) simulation to large systems by reducing the computational effort. The progress in 
these two approaches will enable us to bridge micro- and macro-scale phenomena. In this 
paper we address the latter approach. 

Atomistic simulations should be founded on a quantum mechanical model in order to 
simulate a wide range of materials within a single framework, since the electronic structure 
determines the energy and the forces on atoms. The local-density approximation (LDA) 
to density functional theoryi^i and semiempirical methods such as the tight-binding (TB) 
approximationi'i reduce the complicated quantum many body interaction in condensed mat- 
ter to a single electron problem. The resultant theory has been applied to a variety of prob- 
lems in materials science. However, it exceeds the capacity of modern computers to treat 
large systems that include thousands of atoms, using widely known methods for solving 
the single electron problem such as the conjugate-gradient method, since the computational 
effort scales as the third power of the system size. 

Therefore, several efficient methods with linear scaling algorithms have been proposed 
during the last decade.iHH§ These 0{N) methods can be roughly divided into two categories: 
variational methods and moments-based methods. The former are the density matrix (DM) 
methodSi and the localized orbital (LO) method^i^ which lead to linear scaling algo- 
rithms from the locality of the density matrix. The latter include the bond-order potential 
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(BOP) methodiHll and the Fermi operator expansion (FOE) methodlli0 which are intrin- 
sically linear in the scaling of the computational effort because the enegy and forces are 
expanded in a finite moment expansion. Several applications of these variational methods 
have already been performed for large systems, which have shown the power of these 0{N) 
methodsSii However, several problems remain in these 0{N) methods. 

Firstly, it is well known that the variational 0{N) methods produce large errors in the 
energy of metallic systems with these long range correlations in the density matrix.ii In 
these cases there is no justification for cutting the matrix elements off at short distances in 
the density matrix. Unfortunately, if the cutoff distance is increased to decrease the error 
in the energy, then the calculation effort increases significantly. 

Secondly, it is well documented that within moments-based methods, the vacancy in 
diamond or silicon can not be described within a low number of moments (about 20). 00 
A very large number of moments (about 200) is needed to reproduce the correct vacancy 
formation energy.ii In the BOP method the forces become exact as the bond orders converge 
to the exact values. This implies that the forces are not consistent with the total energy 
if the recursion is terminated at a finite number of levels. In the other moments-based 
methods, such as the FOE methodlilJii and the global density of states method,0S the 
exact forces can be calculated. However, these methods are also unable to reproduce the 
vacancy formation energy within a low number of moments.ii 

Any robust 0{N) method should satisfy the following criteria. Firstly, the method should 
give accurate energies for a wide range of materials (insulators, semiconductors, metals, and 
molecules) with minimum computational effort. Secondly, the Hellmann-Feynman forces 
should be consistent with the total energy at any useful level of approximation. Thirdly, the 
algorithm should be suitable for parallel computation. 

Our goal is to establish the bond-order potential (BOP) method as an 0{N) method that 
satisfies these three criteria. In Sees. 11 and III, we present the theory of the block BOP0 
within the orthogonal tight-binding (TB) representation. We stress that the introduction of 
block elements into the BOP formalism improves remarkably the accuracy of the energy and 



forces. In Sec. IV we explain the reason why the block BOP gives accurate energies with 
a low number of moments by analyzing the vacancy formation energy of diamond carbon 
in terms of the bond-order. In the remainder of this paper the deformation of a single- 
wall carbon nanotube is used to demonstrate the applicability of the method to large scale 
atomistic simulations. 

2. THEORY 

A. Tight-binding 

We develop the block BOP within the two center orthogonal TB represent at ion. @'0 It will 
be assumed that the basis set is an orthonormal set of atomiclike orbitals \ia) where i is 
a site index, and a an orbital index. The Hamiltonian can be represented by the matrix 
Hia,ji3 = {i<^\H\j/3). The on-site elements of the matrix are written as eia- The cohesive 
energy, assuming that the electrons are at a finite temperature T, is the sum of bond, 
promotion, and repulsive energies: 

Ecoh -^bond ~l~ -^prom ~l~ -^repi (l) 

where the repulsive energy is given by the sum of pair potentials or embedded potentials 
which are usually determined so that the TB model reproduces equilibrium structures and 
elastic constants. The bond energy is the attractive contribution that leads to cohesion. 
There are two different but equivalent expressions that describe the bond energy. The first 
gives the bond energy in terms of the on-site density of states as follows: 



where nia{E) is the density of states projected onto orbital \ia), and the function f{x) = 
1/[1 + exp(a;)] is the Fermi function. The second gives the bond energy explicitly in terms 
of the individual intersite bond energies as follows: 




(2) 




(3) 
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where Qia,ji3 is the bond-order between orbitals \ia) and and the expression parenthesis 
represents the corresponding bond energy associated with orbitals \ia) and \ jl3). This allows 
us to interpret the bonding and structure of molecules and solids from a chemical point of 
view.ii It should be noted that the bond-order is not pairwise but is determined by the 
particular arrangement and connectivity of the atoms adjacent to the two atoms forming 
the bond. In the block BOP representation the two different expressions Eqs. (2) and (3) 
for the bond energy are exactly identical at any level of approximation. The proof will be 
given later subsection. The promotion energy is defined by 

i?prom = E (e-^- - ^iNl) ' (4) 
ia 

where Ni^ and A^^^ are the number of electrons in \ia) in the condensed and free atomic 
systems, respectively. The promotion energy is repulsive due to the excitation of electrons 
from their free atomic ground state as the atoms are brought together. Therefore, the 
cohesive energy of a system is determined by the balance between the attractive bond energy 
and the repulsive pairwise/embedding and promotion energies. The bond and promotion 
energies can be repartitioned into the band and atomic energies: 

ia^jP ia 
= -2'band ~ -^atoms- (5) 

-Eband is equal to the energy which is defined by integrating J2iaEnia{E) up to the Fermi 
level. 

In the TB model the single particle eigenfunctions are expanded in a basis set that is an 
orthonormal set of real atomiclike orbitals: \ia). 

I0) = E^S^K«)' (6) 

ia 

where the expansion coefficients are defined by Cf^"* = {ia\(f)). C\'j^ is always real because 
of real atomic orbitals and Hamiltonian. Then the bond-orders may be defined in terms of 
the expansion coefficients as follows: 



= 2EC^j^^C7iiV(^^^), (7) 

where the factor 2 accounts for spin degeneracy, e*^*^-* is the eigenvalue corresponding to an 
eigenstate 

The force on atom k is obtained by differentiating Eq. (1) with respect to atomic posi- 
tions: 



The first term of Eq. (8) is identically zero so that 



where the first term of Eq. (9) is the Hellmann-Feynman force. If the bond-orders are 
approximate values, then the sum of the derivatives of the bond-orders with respect to 
atomic positions will not be zero, so that Eq. (8) gives the exact force which is consistent 
with the total energy. However, in the block BOP representation the forces are given by 
Eq. (9), since it is very difficult to evaluate the derivatives of the bond-orders. Hence, the 
forces calculated by block BOP become exact as the bond-orders converge to the exact 
values. In Sec. 2 the compatibility between the force and the energy will be discussed from 
numerical tests using constant energy molecular dynamics simulations. 

Although the forces are not consistent with the total energy in the usual BOP methods, 
it is possible to evaluate the exact forces at any level of approximate by the other moments- 
based method, global density of states method.^ However the use of the global moments, 
which are introduced to decrease the computational effort, leads to a reduced rate of con- 
vergence of the energy as a function of the number of moments. In Appendix of this paper 
we present a novel method to evaluate the exact forces. 
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B. Block bond-order potential 



The local density of states and bond-orders can be related to the one particle Green's 
functions. The one particle Green's function operator is defined by 

G{Z) = (Z - H)-^ 

Then the imaginary part of the diagonal elements of the Green's function matrix give the 
local density of states: 

Tmr (F I - V -o^(^Q^I'^)('^Kq^) 

Im G^aM^^ + lU J - 1. ^E-eWy + {0+Y 
= -nUiaiE). 

Therefore 



TiiaiE) = --Im GiaM^ + iO+), (11) 

TT 

where Gia,ia{Z) = {ia\G{Z)\ia) , 0^ represents a positive infinitesimal, and 6{x) is the delta 
function. The imaginary part of the off-diagonal elements of the Green's function matrix 
has the following relation to the expansion coefficients of the single particle eigenf unctions: 

Im G,^,p{E + iO+) = -TT 5: C^fC^t^SiE - e^t>)), (12) 

Multiplying the both sides of Eq. (12) by the Fermi function, integrating with respect to the 
energy we obtain the following useful expression for the bond-order: 



lmlG,^,jp(E + iO+)f{^^)dE 

- -^EC^fcli^ I S{E - e^'^)f{^)dE 



.(0) 



Therefore 

Q^aJp = --lm f G,^,,p{E + iO+) f{^-^)dE. (13) 

The evaluations of the bond energy Eqs. (2) and (3) require calculating the local density 
of states and bond-orders. We obtain the local density of states and bond-orders from the 
Green's function through Eqs. (11) and (13). The diagonal elements of the Green's function 
matrix can be calculated in a numerically stable way by the recursion method.@'0 Block 
BOP is a general recursion method for evaluating efficiently both the diagonal and off- 
diagonal elements of the Green's function matrix by the recursion method. The first step of 
the recursion method is to tridiagonalize the Hamiltonian using the Lanczos algorithm. § In 
the block BOP we introduce the block Lanczos algorithm with the starting state as a single 
site containing all the valence orbitals rather than the usual scalar Lanczos algorithm with 
a single starting orbital.ciJ However, the application of the conventional block alg orithmBU 
to finite systems such as molecules introduces a numerical instability, since the terminal 
number of recursion levels of the vr bond are different from that of the a bond in the 
recursive algorithm. Therefore, we modify the conventional block Lanczos algorithm. A 
series of procedures for the modified block Lanczos algorithm can be carried out as follows: 

\Uo) = {\tl),\^2),...,m,)). (14) 

A„ = {Un\H\Un). (15) 

|rn) = H\Un) - |f/„-i)*5„ - \Un)A^. (16) 
{B,,+i? = {rn\rn). (17) 

(An)' = Vn(^n+l)V,. (18) 
5„,+i=An*i:n- (19) 
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(20) 




rn)iRn+l) 



-1 



(21) 



and 5^ are recursion block coefficients with Mj x Mj in size, where Mj is the number 
of atomic orbitals on the starting atom i, and the underhne indicates that the element is a 
block. 

The states |t/„) = (|L„i), |L„2), • " " > l-^nMj)) represent the Lanczos basis, and are or- 
thonormal and block-tridiagonalize the Hamiltonian. The modified algorithm gives different 
expressions for the block elements B_n+i ^^'^ these inverses compared with the conventional 
algorithm. The block elements in the conventional block Lanczos algorithm are defined by 



The failure in the conventional algorithm can be illustrated by a carbon trimer with a 
linear chain structure along the x-axis. If the block Lanczos algorithm is applied with 
the central atom in the trimer as the starting state, then the Py and orbitals span two 
independent subspaces. Thus, the recursive algorithm finishes after only one iteration for 
the Lanczos vectors concerned with the Py and pz orbitals. This gives two zero eigenvalues 
in the four eigenvalues of the block element (^2)^. Then one can not evaluate the inverse 
of B_2 using Eq. (23). Therefore, defining B_2 and its inverse by the modified Eqs. (20) and 
(21), respectively, and assuming that the diagonal elements of A^^ corresponding to the zero 
eigenvalues are zero we have 



Rn+l — Y-n^n Y-n- 



(22) 




(23) 




-1 



1 



^2(^2) 



(24) 







V 
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\ 



1 



(t/2|C/2) 



(25) 







V 







IC/2) is reduced to the state with two vectors, while the starting state \Uo) is constructed by 
the four vectors, which permits us to iterate once more with the recursive algorithm. The 
conventional block Lanczos algorithm does not satisfy both Eqs. (24) and (25), since the 
block elements and the inverse are obtained from the unitary transformations of and 
the inverse, respectively. Therefore, the conventional algorithm terminates at this recursion 
level even though the Lanczos vectors for the a orbital can still hop. This reduction of the 
state avoids the numerically instabilities for the case of small eigenvalues of {B_^^iy, even 
when the eigenvalues are not zero. 

Application of the block Lanczos algorithm defines an orthonormal basis set called the 
Lanczos vector or basis. The Lanczos vectors reflect the neighboring atomic arrangement of 
the starting site. In Fig.l we show the Lanczos vectors on an s-valent square lattice. The 
Lanczos vectors spread gradually from the central atom as the number of recursion levels 
increases. Thus, we now expand a one electron eigenstate using the Lanczos vectors: 



where D^^^ = {Lnu\4')- Then the representation based on the atomic basis can be transformed 
into that of the Lanczos basis set by the matrix U such that 



where U is defined by {ia\Lnu), and T can be the Hamiltonian H, the derivative of Hamil- 
tonian with respect to atomic position dH/dVi, the bond-order 0, or the Green's function 
G{Z) matrix. The index L indicates the representation based on the Lanczos basis. Equa- 
tion (27) is a pseude unitary transformation, and the matrix U becomes unitary when the 
number of the recursion levels is infinity in infinite systems. If the block Lanczos algorithm 



(26) 



UTU, 



(27) 
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is started through Eq. (14) with the atomic orbitals on atom i as the starting state, then 
considering Eq. (27) and the orthonormahty of the Lanczos basis, we can relate the bond 
orders in the Lanczos basis representation to the bond-orders based on the atomic basis by 
the following simple relation: 



(28) 



where Q^j and Go„ are the block elements of the bond-orders for the atoms i and j, and the 
states \Uo) and \Un), respectively. For example 6j,- signifies 



0il,7l 0il,72 



i2,jMj 



(29) 



where Mj and Mj are the numbers of atomic orbitals including atoms i and j, respectively. 
In Eq. (28) *£/„,-, which is the (n,j) block element of the matrix *f/, is defined by 



(L„,i|jl) (L„i|j2) 



{Lnl\jM,) 
{Ln2\jM,) 



\ 



(30) 



The simple relation Eq. (28) allows us to evaluate the bond-order in terms of the Lanczos 
basis representation. We have only to calculate the 0th block line, which are the bond-orders 
between the starting atom and the Lanczos vectors surrounding the atom, of the bond-order 
matrix. In the block BOP the bond-orders are evaluated in the Lanczos basis representation, 
and then we get the bond-orders based on the atomic basis from Eq. (28). 

It is essential to start the block Lanczos algorithm with a single site as in Eq. (14). 
Although it is possible to derive an analogous transformation to Eq. (28) using the usual 
scalar Lanczos algorithm, the bond energy of the system depends on the rotation of the 
system.ii Thus, the use of the scalar algorithm is not appropriate, since the bond energy 
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should be invariant to the rotation of the system. We could also start the recursion with a 
cluster containing a neighbor shell of atoms instead of a single site.H However, this choice 
is unsuitable because it is highly computationally intensive. 

In the Lanczos representation the Hamiltonian is block-tridiagonalized: 



{U^\H\Un 



if m = n. 



*i?„ if m = — 1, 

(31) 

B_n^i if m = + 1, 
otherwise. 



The block element G_qq{Z) = {Uo\G\Uo) can be written explicitly by the form of the multiple 
inverse, since the Green's function matrix G{Z) is the inverse of the matrix {Zl — H). 
Appling repeatedly the partitioning method, which is a method for calculating the inverse 
of matrices, to the matrix {Zl — H) we get 

G^,{Z) = [Zl -A,- 'B,[Zl -A,- 'B^[ ■ r'm-'B.rK (32) 

Gqq{Z) is equal to the block element G_^^{Z) based on the atomic basis, since we have started 
the block Lanczos algorithm with Eq. (14). Therefore, the local density of states can be 
evaluated from the diagonal elements by Eq. (11). Also the trace of Gqq{Z) gives the local 
density of states on atom i. 

Moreover, by taking account of the block-tridiagonalized Hamiltonian and the identity 
{Zl — H)G{Z) = I in the Lanczos basis representation, the off-diagonal elements of the 
Green's function matrix Gg„ may be obtained from the following recurrence relation: 

Gi;SZ) = (GoVi(^)(^I - An-i) - Gt-2iZ)%.-i - 5inl) (^J^S (33) 

where 6 is the Kronecker's delta, Go_i(2') and *5g are 0, respectively. All the off-diagonal 
block elements Go„(Z) are related to the diagonal block element Gqq{Z). Once Gqq(Z) has 
been obtained, the off-diagonal block elements are easily evaluated from the above recur- 
sive relation. The simplicity of evaluating the off-diagonal block elements is an important 
advantage of the Lanczos basis representation. The block elements of the Green's function 

12 



matrix have the same relation to the bond-orders based on the Lanczos basis as that of the 
atomic basis representation: 



= — Im / Go^{E + iO+) f{-—^)dE (34) 

In case the bond-orders are evaluated by Eqs. (28) and (34), we can prove that the 
two different expressions Eqs. (2) and (3) for the bond energy are identical at any level of 
approximation. Consider the trace of G{Z){Z\ — H). Transforming the trace of the atomic 
basis representation into that of the Lanczos basis using Eq. (27), and making use of the 
identity G{Z){Zl — H) = I in the Lanczos basis representation we see that the trace is a 
constant: 

tr {G{Z){Zl - H)} 

= J2^T {ZG,,{Z)} -J2^t{G,^{Z)H^,} 

i ij 

= Etr{^G^cS^'(^)} -Etr{G^r(mn"i'^} 

i in 

= Etr(I«), (35) 

i 

where Ij is a unit matrix with Mj x Mi in size. The index L^*-* indicates the representation 
based on the Lanczos basis with the starting state on atom i. Considering the imaginary 
parts of the trace we have 

lmJ2ZGia,ia{Z) = Im E Gia,jf3{Z)Hjp^ia. (36) 

We see that the two expression for the bond energy give the same energy, since the Green's 
functions can be related to the local density of states and bond-orders through Eqs. (11) and 
(13), respectively. The block BOP, thus, provides the equivalence of the two expressions for 
the bond energy in a natural way, whereas in the usual BOP the Green's functions need a 
carefully chosen truncator in order to satisfy the sum rule.i 
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C. Moment description 



The moments of the local density of states allow us to link the behavior of the electronic 
structure to the local topology about the given site. We now discuss the relation 

between the block recursion matrices and the moments of the density of states. From 
Eq. (10) for \Z\ oo, the diagonal element Gqq{Z) can be rewritten as follows: 

(f/o|0)(0|f/o) 



Z - e(*) 



E4ME 



ZP' 



-1 



E 

p=0 



where 



On 



iOO 



(37) 



y ^il ^ip ^ip ' ' ' ^ip ^ip j 



(38) 



,,(P) - ^ M (M) 



(39) 



and /i^^^ is the block element of the pth moment for the atom i, the diagonal elements of 
which give the pth moments of the projected density of states UiaiE). Thus, Eq. (37) is 
the moment expansion of the Green's function Gqq{Z). Also the pth block moment can be 
evaluated explicitly as the expectation value of the pth power of the Hamiltonian in terms 
of the block elements A„, 5„: 



= J2 iUo\H\u^,)iamAH\u^,)---iu^,,AH\Uo) 

The first few block moments are 



(40) 
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S - {Aof + 'B,B,. (41) 

Prom Eq. (40) we see that the pth moment is the sum over all self-returning paths of length 
p. The first moment corresponds to a hop on a single site, the second to nearest neighbors 
and back, and so on. Thus, the atomic connectivity can be related directly to the electronic 
structure through the description of the Green's function by the moments. 

Multiplying both sides of Eq. (37) by {E + 0+)'', and integrating with respect to the 
energy E we get the following relation: 

-Im / E^^G^,{E + 0+)dE = (42) 



TT J-oo ■ ^0 



This relation means that the imaginary part of the moment of the block diagonal element 
in the Green's function matrix is equal to the moment of the Hamiltonian. 
Let us define the orthogonal block polynomials P„(x): 

xRnix) = Pnix)A^ + En-lixYB^ + Pn+l{x)B^+„ (43) 

where P-i{x) and Po(^) the zero matrix and the the unit matrix I with x Mi 
in size. By using the block polynomials the recursion block elements and 5„ can be 
expanded with the moments: 

A„ = (C/„|i/|C/„) 

= 'P^{H)iUo\H\Uo)P^{H) 

2n+l 

= E Q.ml£^^gLm- (44) 



= {Un\H\Un-l) 

= 'Pn{H){Uo\H\Uo)Pn-l{H) 



2n 



-E^m/fi^^m- (45) 
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In the derivations of Eqs. (44) and (45) we have assumed the substitution: \Uo)H — > H\Uo) 
and H{Uo\ iUo\H. The block coefficients a^, a^, 6^, and 6^ are given by the recursion 
block elements. For example and B^^ can be written as follows: 

4i = - ^HS -/fS'4„ + 4„rfo'A} (46) 

a = ('&)-' {a£'-4„Mii'}- (47) 

In case the recursion in the block Lanczos algorithm is terminated at the gth level, the 
diagonal block element of the Green's function matrix can be expanded with the {2q + l)th 
moments, because it is constructed by the multiple inverse with the recursion block elements 
A^{n — q), B_^{n — 1 ^ q) given by the qth recursion. As shown in Eqs. (44) and (45), 
the recursion block elements are expanded in terms of the moments. Thus, Gqq contains the 
~ (2g + l)th moments. This implies that up to {2q + l)th moment is included in the sum 
of the moment expansion Eq. (37), and Eq. (42) satisfies for r < 2g + 1. 

To obtain the moments for the off-diagonal elements of the Green's function matrix, 
multiplying both sides in Eq. (33) by {E + 0+)'' and integrating with respect to the energy 

we have 

E^G^{E + Q+)dE = E / E^+^G^iE + 0^)dE) c^, (48) 

where the block coefficients can be written in terms of the recursion block elements. As 
mentioned above the right side of Eq. (48) is equal to the moment of the Hamiltonian for 
r + m < 2gf + 1, so that the left side gives the exact moment //^^^ for r < 2gf + 1 — n. This 
means that the off-diagonal elements of the Green's function matrix can be expanded with 
up to the (2g + 1 — n)th moment, which results in the expansion of the bond-order 9q„ by 
up to the (2g -|- 1 — n)th moment. Moreover we can relate the bond-orders in the atomic 
basis representation to the moments through the transformation Eq. (28). In the right side 
of Eq. (28) the bond-order ©q^ for n = q determines the maximum order of the moments for 
the bond-orders based on the atomic basis. So we see that the bond-orders in the atomic 
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basis representation can be expanded with the moments for r < + 1. Thus, in the block 
BOP the off-diagonal elements of the Green's function matrix can be constructed with the 
moments for r < q + 1, while the diagonal elements have the information of the moments 
for r < 2q + 1. This could be imply the difference in the convergence properties of the 
bond energy and the forces. On a simple consideration it is estimated that the rate of the 
convergence of the force is about half as fast as that of the bond energy in terms of recursion 
levels. However it should be noted that the contribution of ©q„ to Q^j decreases as the 
recursion level n increases, since the Lanczos vectors, which hop repeatedly in the atomic 
connectivity, have their weight away from the starting atom as the recursion level n increases. 
Thus, the bond-orders in the atomic basis representation do not have all the moments of 
the higher order more than the (g -|- l)th, but can include the higher moments through the 
for n < q. In this case whereas the inexact moments for r < 2g -h 1 — n are included 
in the bond-order in the atomic basis representation, the error can be negligible, since the 
bond-orders Gq^ become small as the recursion level n increases. So it is stressed that the 
higher moments can be included in the bond-order based on the atomic basis through the 
Green's function Gq„ for small recursion levels n. Therefore, it is expected that the forces 
should be comparable to the bond energy in terms of the convergence rate. In Sec. 3 we will 
discuss this point again numerically. 

D. Details on implementation 

The technical details to implement the block BOP are given in this subsection. For an 
infinite system, there could be an infinite number of levels in the multiple inverse of the 
diagonal Green's function. It is often the case, however, that the exact values can be 
replaced by estimated values after a certain number of levels, without reducing the accuracy 
significantly. The simplest approximation is to take A^^ = A^, 5„ = for n > rit, 
where rit is the number of exact levels, and and are constant block elements. This 
approximation is reasonable from the observation that the scalar elements in both A„ and 
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B_„ converge to constant values or oscillate around constant values as n tends to infinity.ii 
We have only to replace the level for n = rit + 1 in the multiple inverse with the terminator, 
since the constant terms can be summed exactly. The terminator can be written by a closed 
form including itself as follows: 

T(Z) = [Zl -A^- 'B^T{Z)B^]-\ (49) 

However, this is still a difficult set of equations to solve, so to simplify matters we assume 
that the off-diagonal elements of T(Z) are zero and all the diagonal elements are the same, 
since the differences between the diagonal elements of and 5„ become small as the 
number of the recursion levels increases, respectively. Then the identical diagonal element 
t{Z) of T_{Z) is written as the square root terminator: 



t{Z) = [Z-a- hH{Z)]-^ 



(50) 



2b V V 26 

where a and 6^ are given by the means of the diagonal elements of A^^ and B_^^ , respectively. 
Thus, we see that the effect of the terminator is to smear out the sharp states with energy 
a into semielliptical bands. The degree of smearing is given by b. 

There are two ways to conserve charge neutrality in the system: local charge neutrality 
(LCN)@ or the total charge neutrality (TCN).t3 Within LCN the on-site energies are varied 
(keeping the splitting between on-site s and p energy levels fixed) in order to conserve the 
number of electrons on each atom. If the excess charge on site i is Qi = Zi — J2a ^ia, where 
Zi is the effective core charge, then the on-site energies can be shifted using the response 
function Xi = J2a -^ia for atom i as follows: 

= ^ia - A^, (51) 

where A is a parameter to accelerate the convergence, and generally is 1.0. The response 
function projected on an atomic orbital ia is given by 



= -Im / [G,^,UE + tO+)ffi^—^)dE. (52) 



Usually no more than three or four iterations are required to achieve the convergence so 
that the absolute value of Q/atom is below 10~^, since — dNia/deia- The assumption of 
LCN has the advantage that the Madelung energy contribution is zero, so that the TB model 
needs not take this into account in its expression for the energy. Also LCN is suitable for 
parallel computation, since the calculations of the bond energy and the forces of each atom 
are perfectly independent within the assumption. However, LCN brings an inefficiency in 
terms of computational effort, since LCN requires the Lanczos algorithm to be implemented 
again, after the charge neutralities of all the atoms has been achieved, since the recursion 
block elements are varied by the shift of the on-site energies. Thus, the block Lanczos 
algorithm and the shift of the on-site energies must be repeated until self-consistency is 
accomplished. This self-consistency requires typically twenty iterations. This discourages 
us from applying LCN in the molecular dynamics simulations. On the other hand, we can 
conserve the total number of electrons in the system by a shift of the chemical potential in 
terms of TCN. If the excess charge of the system is Q — Qh then a good approximation 
of the chemical potential is given by 

/.' = /. + a|, (53) 

where X = J2i^i- The convergence is achieved after only three or four iterations. The TCN 
assumption, corresponding to the micro canonical distribution, has physically appropriate 
meaning, which is consistent with the usual electronic structure calculations by diagonaliza- 
tion. Moreover within TCN we need not repeat the Lanczos algorithm, since the recursion 
block elements are not varied by the shift of the chemical potential. Thus, TCN has consid- 
erable advantage in terms of computational effort. The TCN condition reduces the separa- 
bility of individual atoms in the calculations of the band energy and forces, and comphcates 
slightly the parallelizability of the program code. However, the evaluation and integration of 
the Green's function, which are time-consuming steps, are performed separately. Therefore, 
we use the TCN constraint to conserve the total number of electrons. 

It is required to integrate the Green's functions with the Fermi function in order to eval- 
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uate the bond energy, bond-orders, and response functions. The integration can be carried 
out in the complex plane by summing up an infinite series over the modified Matsubara 
poles.0'0S The general form can be given as follows: 



where A{x) is an arbitrary function defined in the complex plane, and (3 = l/ksT. Also 
Ep are the poles of the approximated Fermi function in the complex plane. This modified 
Matsubara summation converges rapidly with about 40 complex poles (P ~ 40) with a 
high electron temperature {ksT > 0.1 eV), although many poles are needed to achieve the 
convergence with a lower electron temperature. In the case of systems with a gap between 
the valence and conduction bands, we need to pay attention to the evaluation of the chemical 
potential, since the response functions in the gap become zero as ksT tends to 0, so that it is 
difficult to estimate the chemical potential under a low electron temperature using Eq. (53). 
This can be solved by smearing the density of states under a high electron temperature. 
Thus, it is required to evaluate the response functions at high electronic temperatures in 
order to obtain stable MD simulations. 

We now estimate the time-dependence within the block BOP. The total system is divided 
into finite clusters centered on individual atoms in order to evaluate the energy and force of 
each atom. The size of the finite cluster is not determined by the size of the total system, 
but by the system and the condition of the MD simulation. Therefore, the computational 
effort is proportional to the number of atoms A'atonn so that the number of computational 
operations can be written as cAi'atom, where c is a proportionality constant. The scaling of 
the constant c can be estimated as a function of the numbers of recursion level q, atoms 
within a finite cluster Uc, and orbitals on an atom M. For simplicity it is assumed that the 
system consists of only one type of element with M orbitals. In the block Lanczos algorithm 




p-i 



(54) 



with 




(55) 
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the time-consuming step is the product of the Hamiltonian matrix by the vector, so that the 
count of operations in the block Lanczos algorithm is nearly proportional to qn^M. At the 
next step, the inverses and recursive calculations are required to evaluate the diagonal and 
off-diagonal elements of the Green's function matrix, respectively, and their integrations are 
performed as the sum of the residues for the poles in the complex plane, so that the count 
of operations for the evaluations is almost proportional to qPM^. Thus, the proportionality 
constant c can be estimated as x qnlM + cq x qPM^, where cl and cq are pref actors 
of the count of operations for the block Lanczos algorithm and the the evaluation of the 
bond-orders, respectively The prefactors depend on the computer, and the system, and the 
criterion of charge neutrality. For example, for the case of a 3 hop cluster, 10 recursion 
levels, and 40 complex poles for diamond carbon, the calculation time of the block Lanczos 
algorithm is comparable to that in evaluating and integrating the Green's functions. 

In the remainder of this section the procedure for implementing the block BOP is enu- 
merated. (I). The partition of the system. The hopping range of each atom is determined 
by terminating the system. There are two ways to terminate the system. One of them is 
the physical truncation that the terminated cluster contains atoms within a sphere with a 
certain cutoff radius. The physical truncation can bring inaccurate properties into the con- 
vergence of the energies, since atoms that have no bonding to other atoms can be included 
in the neighborhood of the cluster surface. Moreover, in MD simulations the energies can 
jump discontinuously when an atom moves in or out of the surface of the sphere. The more 
stable way is logical truncation. The cluster of size n is here defined by all neighbors that 
can be reached by n hops. Provided the cutoff distance for the hopping integral is identical 
to that defining the connectivity of the bonding, the energies are continuous as a function 
of time in MD simulations. Therefore, it is desirable to truncate logically the system in 
terms of accuracy. (11). The block Lanczos algorithm. The Hamiltonians for the individual 
terminated clusters are constructed. For these small cluster Hamiltonians the block Lanczos 
algorithm Eqs. (14)~(21) is applied. (III). The evaluations and integrations of the Green's 
functions. In the Lanczos basis representation the diagonal and the off-diagonal elements 
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of the Green's functions are evaluated using Eqs (32) and (33), respectively, and then their 
integrations are performed via the modified Matsubara summation with Eq. (54). (IV). The 
transformation into the atomic basis representation. The bond-orders based on the Lanczos 
basis are transformed into those in the atomic basis representation using Eq. (28). (V). The 
bond energy and forces. From Eqs. (3) and (9) the bond energy and forces are evaluated, 
respectively. 



3. CONVERGENCE PROPERTIES 

0{N) methods with linear scaling algorithms are approximate approaches compared to the 
exact diagonalization for dealing with large scale systems, so that the realization of the 
0{N) algorithms is accompanied by decreases in computational accuracy in exchange for 
computational efficiency. Therefore, 0{N) methods should only be applied to atomistic 
simulations once their accuracy and efficiency has been tested. 

In the block BOP two approximations are introduced to reduce the computational effort: 
the number of moments, or recursion levels, and the size of the cluster of atoms over which 
the hops are made. The finite approximations for the number of levels and the size of the 
cluster can lead to the errors in the energies and forces. Thus, we now investigate the block 
BOP through several test calculations in terms of its accuracy and efficiency. In order to as- 
certain applicable bounds for a wide range of materials, the energy and force convergence are 
examined for an insulator (carbonil in the diamond structure), a semiconductor (siliconll), 
a metal (titanium, described by a canonical d-band model), and a molecule (benzene^) as 
functions of the number of recursion levels and the size of cluster. Moreover, in terms of 
the computational efficiency the block BOP is compared with k-space calculations in com- 
putational time. Also as a test of the quality of the forces, we perform a constant energy 
molecular dynamics (CEMD) simulation of carbon. Figure 2 shows the cohesive energy per 
atom for carbon in the diamond structure, silicon in the diamond structure, hep titanium, 
and benzene. The cohesive energies were calculated using 2 ~ 15 recursion levels for three, 
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five, and seven shell clusters by the logical truncation method, where the three, five, and 
seven shell clusters for the diamond structure include 41, 147, and 363 atoms, respectively, 
and these clusters for the hep structure contain 153, 587, and 1483 atoms, respectively. The 
cohesive energies for carbon and silicon converge rapidly to the results of k-space calcula- 
tions. The errors for carbon and silicon are only 1 % at six recursion levels. Thus, we see 
that up to the 13th moment corresponding to six recursion levels determine the cohesive en- 
ergies. The contribution of the higher order moments is unimportant, since the convergence 
properties are almost identical for three, five, and seven shell clusters. The cohesive energy 
for silicon converges more slowly compared with that of carbon in the rate of convergence 
for the size of cluster. This suggests that a semiconductor such as silicon requires higher 
moment than an insulator such as carbon for good convergence of the cohesive energy. The 
cohesive energy for the metallic hep titanium converges very quickly in terms of the number 
of recursion levels. For the five and seven shell clusters the cohesive energy converges fully 
to the k-space result, while the convergence value for the three shell cluster is in error by 2 
% from the k-space result. For benzene the convergence is achieved with a very small cluster 
(2 shells). The error at four recursion levels is only 0.1%. We see that the block BOP can 
evaluate accurately the cohesive energy for a molecule with a sparse structure like benzene, 
which has both localized a bonds and delocalized vr bonds. 

The calculation of the vacancy formation energy is a severe test to distinguish the ac- 
curacy of different 0{N) methods, since it is a criterion that tests the precision which the 
dangling bonds caused by the vacancy are handled by 0{N) method. In practice, the usual 
moment-based 0{N) methods fail to reproduce the vacancy formation energy of carbon in 
the diamond structure even when dozens of moments are included.ii@ The computational 
error at 30 moments is still about 20 % compared to the k-space result. In Fig. 3 we show 
the vacancy formation energy for carbon in the diamond structure, silicon in the diamond 
structure, and hep titanium. These are calculated as the difference between the energy for a 
bulk unit cell (of 64, 64 or 32 atoms, respectively) with a single atom removed, and the per- 
fect bulk cell energy scaled to 63, or 31 atoms. The results are for an unrelaxed vacancy. The 
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convergence properties for carbon and silicon are almost identical. The vacancy formation 
energy in the five and seven shell clusters converges smoothly toward the k-space results, 
while in the 3 shell cluster the converged values for carbon and silicon are 15 %, and 13 % 
underestimated, respectively. In the seven shell cluster at 15 recursion levels the errors for 
carbon and silicon are only 1%. Thus, we see that the block BOP gives an accurate vacancy 
formation energy for strongly covalent materials such as carbon and silicon with the use of 
about 30 block moments. This remarkable result suggests that the block BOP accurately 
describes dangling bonds in comparison with the usual moment-based methods. In Sec. 4 we 
will discuss the advantages inherent in the block BOP by analyzing the vacancy formation 
energy in terms of different bond-order contributions. For titanium the vacancy formation 
energy converges to the k-space result equally within the three, five, and seven shell clusters. 
The error for the 3 shell cluster at 5 recursion levels is about 6%. The vacancy formation 
energy oscillates with respect to the number of recursion levels due to the long range value 
of the density matrix (see fig. 2 of ref. 23). The oscillations are damped by imposing LCN 
instead of TON to conserve the number of electrons. 

The accuracy of the forces is investigated from two different perspectives. The first 
is the accuracy when compared to the exact k-space result, the second is the degree of 
correspondence between the numerical and analytic Hellmann-Feynman forces. In order to 
perform rehable MD simulations the two criteria should be satisfied. In Fig. 4 we show 
the z-component of the force on an atom in the bulk-terminated (001) surface of carbon, 
silicon, and hep titanium, and the force on a hydrogen atom on benzene. For carbon the 
force of the three shell cluster overestimates by about 130 % in comparison with the k- 
space result, although the error in the Hellmann-Feynman term is only 1 %. The forces of 
the five and seven shell clusters converge smoothly toward the k-space result. The rate of 
convergence in silicon is much better than that of carbon. Even the three shell cluster shows 
a converged value that differs by only 5 % from the k-space result. The three, five, and 
seven shell clusters of Ti show similar convergence properties of the forces, the converged 
value being underestimated by about 8% compared with the k-space result. For benzene the 
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force converges rapidly with small cluster size. As discussed in Sec. 2 the bond-orders can be 
expanded using the lower order moments compared with the density of states in the block 
BOP. It can be estimated that the forces should converge more slowly at the k-space results 
than the bond energies, since the forces on the atoms are evaluated using the bond-orders. 
However, these numerical results for the forces show that the convergence rate of the force 
is comparable to that of the bond energy. This means that the sum of Eq. (28) converges 
rapidly as the number of the recursion levels increases because of the diffusion of the Lanczos 
vectors. 

As a test of the consistency between the total energy and the forces, constant energy 
molecular dynamics (CEMD) simulations have been performed for carbon. If the forces are 
equal to the derivative of the total energy with respect to atomic positions, the total energy 
of the system is conserved. Thus, the CEMD simulation is a criterion to investigate the 
consistency of forces. In Fig. 5 we show the energy for carbon at 1000 and 5000 K as a 
function of time using five and ten recursion levels. The initial structure is the diamond 
lattice, and the unit cell is fixed in volume and shape. When the initial temperature of 
the system is 1000 K, the atoms oscillate around the equilibrium positions. At five and ten 
recursion levels we see that the total energy is almost conserved. When the temperature 
is raised to 5000 K, the carbon in the diamond structure transforms into liquid carbon 
with mainly three coordinate structure. From Fig. 5 we see that the forces are of good 
quality at ten recursion levels, while the total energy at five recursion levels increases by 
about 10 eV during the 1 ps, which corresponds to a temperature increase of 1800 K. These 
results indicate that the block BOP can give forces consistent with the total energy, provided 
the proper number of recursion levels is used, even for liquid materials such as carbon at 
a high temperature. On the other hand, in the variational DM method, although only 
the Hellmann-Feynman term survives formally as the derivatives of the band energy with 
respect to atomic coordinates, total energy of liquid silicon in the CEMD simulation exhibits 
a steady upward drift .S 

To study the computational efficiency of the block BOP we carry out two benchmark 
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tests: the comparison between the block BOP and the k-space calculation in computer time, 
and the relation between the computational error and the computer time. Figure 6 shows the 
time to evaluate the energy and forces for a cell containing carbon in the diamond structure 
as function of the number of atoms in the cell for the block BOP and k-space using a single 
k-point. The crossover point at which the block BOP becomes favorable is about 100 atoms. 

Figures 7(a) and (b) show the relation between the error and the the time per atom 
to evaluate the energy and forces in the calculations of the vacancy formation energy of 
diamond carbon and hep titanium, respectively. We see that the block BOP can calculate 
the vacancy formation energy to high accuracy within almost the same computational time 
as the other moment-based results reported by D.R.Bowler et al.il where the calculations 
were performed using the same computational facilities. We note that the block BOP has 
given a good convergent result of the vacancy formation energy in diamond carbon for 
the first time with a moments-based method, while the computational time to achieve this 
convergence is still ten times slower than that of the DM method. This work, therefore, 
still supports the conclusions of the study in ref. 23 that the DMM is best for systems with 
energy gaps, but that moments-based methods such as BOP are best for metallic systems. 



4. ANALYSIS OF VACANCY FORMATION ENERGY 

The block BOP can provide chemical insight into the nature of the bonding in molecules 
and solids in terms of the bond-order. The bond-order is an useful quantity indicating the 
strength of bonding between two atoms. In practice, it is well known that the bond length 
is nearly proportional to the bond-order for the vr bonded hydrocarbons. 

In this section we analyze the vacancy formation energy of carbon in the diamond struc- 
ture in terms of the bond-order, and discuss the reason why the usual moment-based methods 
can not reproduce the vacancy formation energy even with several dozens of moments.ilifl 
The reduced TB mo is introduced in order to clarify the analysis of the vacancy 

formation energy in terms of the two scalar bond-orders Go- and G^r, respectively. In the 
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reduced TB model the three independent bond integrals hssa, hppa, and hgpa are reduced 
to the two independent variables h„ and p„ by assuming that h^p^^ is the geometric mean 
of \hsscr\ and /ippo-. This allows the a bond energy between atoms i, and j to be written as 
the single quantities h„{Rij)Q'^J rather than the sum of three terms hssa^ssa, '^hspa^spa, and 
hppa-Qppa, respectively. That is, we can write 

-^bond — -^o- + -^TT 

= -2/i^'^) e^'^) - 4/i(;^) ei'^) , (56) 

where 

= {I + Pa)\hsscT\, 



^ "sso- I "\ifa"spa i fa^ppa lr<-,\ 
®^ = ' ^'^^ 

with Per = ^pptr/l/isso-l- All the hopping integrals and bond-orders are defined as positive 
quantities. In addition, the cutoff distance in the hopping integrals and repulsive potential 
is reduced from 2.6 A in the original TB fitii to 2.5 A. This modification simplifies the 
analysis, because atoms on the diamond lattice do not interact with this second neighbors 
who lie at a distance of 2.517 A. Also we apply LCN to the analysis rather than TCN, 
since chemical concepts like the promotion energy require the total number of electrons on a 
given atom to be invariant as the atoms we brought together to form the bond. In table I we 
give the cohesive energy and vacancy formation energy of carbon in the diamond structure 
calculated using the original TB and the reduced TB methods. The changes in the cohesive 
and vacancy formation energies introduced by the reduced TB simplifications are only 0.1 % 
and 3 %, respectively. Therefore, it is an excellent approximation to analyze the vacancy 
formation energy using the reduced TB method with LCN. 

Figure 8 shows the diamond lattice with a vacancy. There are four first (FN) and twelve 
second (SN) neighboring atoms about the vacancy. The twenty four third neighboring atoms 
in total are grouped into two kinds of atoms (TN, TN'), each of them including twelve atoms. 
The number of valence s and p electrons and the corresponding promotion energy of the FN, 
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SN, TN, and TN' atoms are given in table II. The number of valence s electrons on the FN 
atom increases by about 6 % compared with that of a carbon atom in the perfect structure, 
which corresponds to an increase of 0.27 electrons in total over the four FN atoms. This 
increase in the number of s electrons on the FN sites reflects that the s component of the 
dangling bond is attracted firmly at the core of the carbon. The number of s electrons on the 
SN, TN, TN' atoms, on the other hand, is almost the same as that of the perfect structure. 
We see, therefore, from table II that 97 % of the total change in promotion energy resides 
in the FN shell of atoms about the vacancy, so that the redistribution of s and p valence 
electrons occurs mainly within the first shell. The change in the promotion energy stabilizes 
the vacancy by 1.858 eV. 

Table III shows the bond-orders and bond energies for a and tt bonds between the pairs 
of atoms FN-SN, SN-TN, and SN-TN', respectively, in the presence of a vacancy. The a 
bond-order for FN-SN increases by 0.4 %, whereas for SN-TN and SN-TN' bonds it decreases 
by 0.6 % and 0.1 % respectively compared with that of the perfect structure. This oscillatory 
behavior in the variation of the bond-orders reflects the screening of the vacancy. However, 
the very small variation in the a bond-order reflects the localized nature of the a bonding in 
carbon which is a strongly covalent material. For the tt bonding the bond-order between FN 
and SN atoms increases by 36 % compared with that of the perfect structure. This increase 
means that the p electron of the dangling bond participates in the tt bonding between 
FN and SN atoms rather than being attracted solely to the core of the carbon vacancy. 
If we had assumed that the bond-orders are invarint to the formation of a vacancy, then 
the vacancy formation energy would have been overestimated by 1.752 eV. This additional 
stabilization energy to the formation of the vacancy is distributed between the a and tt 
bond energies, as -0.799 and 2.551 cV contributions, respectively. Thus, the absolute ratio 
of the (7 to TT contributions is about 1 to 3, which is considerably larger than the ratio of 
the a and tt bonding energy (18.054:0.635) in the perfect diamond lattice. In the a bond 
energy the contribution of the SN-TN and SN-TN' bonding to this stabilization energy is 
comparable to that of the FN-SN bonding. On the other hand, the stabilization energy for 
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the TT bonding comes from mainly the FN-SN bonding as the p electron in the dangling bond 
only participates in the vr bonding between the FN and SN atoms. Finally, the total vacancy 
formation energy (9.783 eV) can be separated as the difference of the repulsive energy (- 
23.983 eV), the bond energy of the absent bonds reproduced by the vacancy (37.380 eV), the 
stabilization of the promotion energy (-1.858 eV), and the stabilization of the bond energy 
(-1.752 eV). 

Figure 9 shows the errors in the bond-orders for cr and vr bonds between the FN and SN 
atoms. We see that the rate of convergence with respect to the number of recursion levels 
of the vr bond-order is twice as large as that of the a bond-order. 

Thus, we have found that the block BOP can separate the different behavior of a and tt 
orbitals correctly. In particular, it can reproduce the different magnitude of reconstruction 
for the vacancy and the convergence rate with respect to the number of the recursion levels. 
In the usual scalar moment-based methods the a and vr orbitals are not separated, since 
an averaged moment is used for the two kinds of orbitals.i This means that the different 
properties of the s and p electrons in the dangling bond are averaged with respect to the 
vacancy formation energy and the convergence rate. As a result a great many moments are 
required in order to reproduce the vacancy formation energy in the usual scalar moment- 
based method. 

5. LARGE SCALE SIMULATION 

The block BOP is applicable to the atomistic simulations of large systems including thou- 
sands of atoms. In this section we discuss the parallel computation required to perform 
such large scale atomistic simulations and illustrate the method with an application to the 
deformation of a single-wall carbon nanotube under compression. It is very easy to paral- 
lelize the program code because of the highly independent nature of the algorithm which 
evaluates the energy and force for each atom. We have only to parallelize essentially the 
three main loops in the program code: the block Lanczos transformation, the determina- 
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tion of the number of electrons on each atom, and the evaluation of forces. In these loops 
independent calculations can be performed for each atom, since no information needs to be 
passed between the individual atoms. The majority of the computational effort is occupied 
by the calculations in these three loops. Thus, if the three loops are parallelized, we have 
almost the ideal parallel code. Figure 10 shows the time to evaluate the energies and forces 
of a diamond unit cell including 512 atoms as a function of the number of processors. It 
is found that the scaleability is almost ideal. The parallelization was done using an auto- 
matic parallel compiler. The ideal scaleability brought by the use of the automatic parallel 
compiler indicates the simplicity of the algorithms within block BOP. 

We have performed a large scale atomistic simulation for the deformation of a single- 
walled (10,10) nanotube under compression along the shaft as an application of the block 
BOP. The nanotube, which has a length of 140 A, includes 2280 carbon atoms. The compres- 
sion along the shaft was performed by the following geometric optimization with a constraint. 
First, the 2;-coordinate of all the atoms in the nanotube oriented along the z-axis are scaled 
as the length of shaft decreases 0.1 % of its initial length. Second, the scaled structure is 
optimized geometrically with a constraint that the ^-coordinate of the atoms within 7 A of 
both ends are kept fixed. By applying repeatedly the optimization to the nanotube, the 
shaft of the nanotube can be compressed statically. In the early stage of the compression 
the nanotube shrinks, maintaining the shape. However, the nanotube buckles periodically 
when the length of shaft reaches about 80 % of the initial structure. Figure 11 shows a snap- 
shot of the ripple-buckling nanotube. The mean wavelength of the ripple-buckling is 4.8 A. 
The appearance of the ripple-buckling is very similar to the behavior observed by transmis- 
sion electron microscope (TEM) and atomic force microscope (AFM) measurements. B§ A 
detail of discussion of the deformation and elastic properties of carbon nanotubes will be 
presented elsewhere. 
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6. CONCLUSIONS 



We have presented the theory of the block BOP based on the Lanczos basis representation 
and the block Lanczos algorithm with a single site as the starting state within the orthogonal 
tight- binding representation. The efficient 0{N) algorithm provides a general recursion 
method for evaluating the bond energy and forces. In the Lanczos basis representation the 
off-diagonal block elements of Green's function matrix can be related to the diagonal block 
elements through a simple recurrence relation. As a result the bond-orders can be easily 
evaluated. Prom the convergence properties for the bond energies and forces it is found that 
the method is applicable to a wide rage of materials (insulators, semiconductors, metals, 
and molecules) with a considerable reduction in the computational effort compared to k- 
space methods. The algorithm becomes more efficient than the k-space calculation when 
the number of atoms exceeds about 100. Constant energy molecular dynamics simulations 
for carbon show that the forces are consistent with the total energy, even if the method is 
applied to hquids. Moreover, the use of the block Lanczos algorithm guarantees that block 
BOP represents the different properties of the a and tt bonds correctly, so that the vacancy 
formation energy of diamond is reproduced correctly for the first time by a moments-based 
method. Pinally, block BOP is very easy to parallehze, the parallehzation of the three 
main loops giving almost the ideal scaleability. The deformation of carbon nanotube under 
compression was demonstrated as an application of the method to large scale atomistic 
simulations. Thus, we conclude that the block BOP is an efficient 0{N) method to perform 
large scale atomistic simulations of a wide variety of materials. 
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Appendix A: 

We derive a novel method for evaluating the exact force, which is consistent with the total 
energy, at any level of approximation. In this method the derivatives of the diagonal Green's 
functions can be evaluated indirectly by making use of the Lanczos basis representation, 
while the method is similar to the global density of states (GDOS) methodH as regards the 
diagonal Green's functions are differentiated. The contribution from the band energy to the 
force is written using Eqs. (2) and (5) as follows: 

-|7,(band) _ 

2^ ^ idG{E + iO 
= — Im / tr 



TT 



/t..|f!y^|B/WdE, (Ai) 



The trace of the right side in Eq. (Al) can be divided into the diagonal block elements for 
individual atoms: 



I OGjE + iO+) j - ^ tr I ^---^^ + 



where the index L^'^^ represents the representation based on the Lanczos basis with atom i 
as the starting site. In the Laczos basis representation the off-diagonal block elements of 
the Green's function matrix are evaluated through the recurrence relation Eq. (33) which 
is derived from the identity {Zl — H)G{Z) = I. Thus, the following useful relation can be 
derived for the derivative of the Green's function: 

= G^"\Z)^^G^'\Z), (A3) 
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Taking account of the (0,0) block element of both sides in Eq. (A3) we have 

Substituting Eqs. (A2) and (A4) for the trace in Eq. (Al) we can derive exactly the con- 
tribution from the band energy to the force at any level of approximation as the following 
very compact form: 

p(band) ^ — E E tr < r^'*'' 



-nm 




i ma,n(3 



where 



= --Im / EG^\e + iO+)G^^'(£; + iO^) f{x)dE, (A6) 

where the summation for atom i is taken for all the atoms which include the atom k in the 
truncation of system. This expression for forces has a very similar form to the Hellmann- 
Feynman force (see Eq. (9)). However r is a new dimensionless quantity indicating the 
strength of bonding between two Lanczos vectors. 

A simple example is shown. Consider the force on atom 1 being at the end of a linear 
s-valent trimer along the x-axis. It is assumed that the trimer has two electrons and the 
hopping integral between the pairs of atoms 1-2, 2-3, and 1-3 are —hi, — /i2, and zero, 
respectively, where the hi and /i2 are positive and h2 < hi, and also on-site energies for all 
the atoms are zero. If the recursion is approximated at the first levels, then the Green's 
functions are given as follows: 

<'^'^-zh^i^zfhi ^^'^ 

<'^'^-z^i^^i- ^ 

G^^\Z) = ^_ + (A9) 

Z-^hl + hl Z + ^hj + hl 
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Gjr(Z) = f^= + ^^=- (AlO) 

Z-^hl + hi Z+^hl + hl 

The band energy of the trimer can be evaluated from the residues of the poles, which are 
occupied, in the diagonal elements of the Green's functions, namely 



£^band = -hi - yjhl + hi (A13) 

In the force evaluation by Eq. (A5) we have only to calculate the TqI^^ (or Tq^^), and TqI^^ (or 
TiQ^^), since the derivatives of Hamiltonian with respect to the position of atom 1 corre- 
sponding to other r's are zero. The integrands for these r's are 

1 -1 i/i 

_ 4 I ~ I 4^1 , 4^1 



Z-hi Z + hi {Z-hif {Z + hif 



ZG^^\z)G^r {Z) 

1 -1 

_ 4 _|_ 4 



(A14) 



z-^hi + hi z + ^hi + h 



+ , + , (A16) 



(z-Tftf + ftl) (z + ^ftf+Ai 



Prom the residues of the poles which are occupied we see that Tqi^\= t[i^^^), and T(f/^^(= 
Tqi^^) are —1/2 and —1/2, respectively. Transforming the derivative of the Hamiltonian by 
Eq. (27) into the the Laczos basis representation we have 

(A16) 
(A17) 

Thus, we can evaluate the contribution from the band energy to the force on atom 1 using 
Eq. (A5), namely 
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dhl 


dxi 


dxi 


dxi ' 

hi dhl 


dxi 


dxi 


yjhl + hi 



ZTi(band) _ o^Lm dH^Q ^ _ r (2) dH^Q ^ 

~ dxr dx, 

-1 dhi -1 hi dhi 

— —2 X — X 2 X — X . - — 

= + (A18) 

The result is identical to the contribution from the derivative of Eq. (A13) with respect to 
the x-coordinate of atom 1. We see that the force by Eq. (A5) is certainly exact. 
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FIGURES 

Fig. 1. The Lanczos vectors on the s-valent square lattice, a), b), and c) are an initial state 
|Lo), \Li), and \L2), respectively. The diameter of the circles is proportional to the magnitude of 
the expansion coefficient in the Lanczos vector. 

Fig. 2. The cohesive energy for carbon in the diamond structure, silicon in the diamond 
structure, hep titanium, and benzene as a function of number of recursion levels for three, five, 
and seven shell clusters, calculated using a square root terminator, and ksT = 0.1 eV. 

Fig. 3. The vacancy formation energy for carbon in the diamond structure, silicon in the 
diamond structure, and hep titanium for three, five, and seven shell clusters as a function of 
number of recursion levels, calculated using a square root terminator, a total charge neutrality, 
and ksT = 0.1 eV. 

Fig. 4. The z-component of the force on an atom on the carbon (001) surface, silicon (001) 
surface, titanium (001) surface, and on a hydrogen atom in benzene for three, five, and seven shell 
clusters as a function of number of recursion levels, calculated using a square root terminator, total 
charge neutrality, and ksT = 0.1 eV. 

Fig. 5. The potential, kinetic, and total energies as a function of time for molecular dynamics 
simulations of carbon using a three hop logically truncated cluster, a square root terminator, total 
charge neutrality, and ksT = 0.1 eV. In panels (a) and (b) the results are for five and ten recursion 
levels at 1000 K, respectively, whereas in panels (c) and (d) they are for five and ten recursion 
levels at 5000 K, respectively. The time step is 0.5 fs. 

Fig. 6. The time to perform the energy and the force evaluation for carbon in the diamond 
structure as a function of number of atoms in the cell for the block BOP, calculated using a three 
hop logically truncated cluster, and k-space. The calculations were performed on an IBM RS/6000 
workstation. 
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Fig. 7. The error in the carbon (a) and titanium (b) vacancy formation energies against the time 
taken per MD step per atom for three, five, and seven shell clusters. The calculations were carried 
out with a square root terminator, total charge neutrality, and ksT = 0.1 eV on a HP9000/735 
workstation. 

Fig. 8. Diamond lattice with a vacancy. FN and SN label the first and second neighboring 
atoms to the vacancy, respectively. Two kinds of third neighboring atoms are distinguished by TN 
and TN'. 

Fig. 9. The errors in the bond-order for a and tt bonds between the FN and SN atoms, where 
the 15 recursion level bond-orders O^- = 0.9156 and Gjr = 0.1412 were taken as the exact values. 

Fig. 10. The calculation time to evaluate the energy and the forces of a 512 atom carbon cell 
as a function of the number of processors by the parallel code. The benchmarks were performed 
on a Sun Ultra 10000 StarFire, using a three shell cluster, five recursion levels, and a square root 
terminator. 

Fig. 11. Ripple-buckling single-wall (10,10) nanotube under compression along the shaft. The 
snapshot is at 80 % of the initial length (140 A). The calculation was performed with a five hop 
logically truncated cluster, ten recursion levels, a square root terminator, total charge neutrality, 
and ksT = 0.1 eV. Within the calculation conditions the error is about 1 % in the bond energy 
for the initial structure compared with the k-space result. 
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TABLES 

Table 1. Comparison of the original and the reduced TB method with respect to the predicted 
cohesive energies of carbon in the perfect diamond structure and tlie diamond unit cell including 
64 atoms with a single atom removed. The calculations were performed with a logical truncation 
of a 7 shell cluster, 15 recursion levels, a square root terminator, local charge neutrality, and 
UbT = 0.1 eV. 





Perfect 


Vacancy 


Vacancy Formation 




(eV/atom) 


(eV/atom) 


Energy (eV) 


Original 








K-Space 


-7.251 


-7.091 


10.110 


BOP 


-7.249 


-7.090 


10.004 


Reduced 








K-Space 


-7.254 


-7.098 


9.860 


BOP 


-7.256 


-7.100 


9.783 
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Table 2. The number of valence s and p electrons and the promotion energies of the first (FN), 
second (SN), and third (TN, TN') neighboring carbon atoms about a vacancy in the diamond 
structure. AE^prom is defined as the difference between the structure with a vacancy and the 
perfect structure in the promotion energy. 





Total 


Perfect 


FN 


SN 


TN 


TN 


Others 


Ns 




1.203 


1.271 


1.203 


1.202 


1.205 




Np 




2.797 


2.729 


2.797 


2.798 


2.795 




-Eprom(eV/atom) 




5.334 


4.886 


5.342 


5.344 


5.328 




A-Eprom (eV/ atom) 






-0.448 


0.008 


0.010 


-0.006 




Number of atoms 


63 




4 


12 


12 


12 


23 


A£;prom(eV) 


-1.858 




-1.793 


0.097 


0.115 


-0.076 


-0.201 


Ai^tSlTotal) X 


100 




96.5 


-5.2 


-6.2 


4.1 


10.8 
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Table 3. The bond-orders and the bond energies for a and tt bonds between the pairs of atoms 
FN-SN, SN-TN, and SN-TN', respectively, in the presense of a vacancy. In the calculations of the 
bond energies and E^^, the hopping integrals are = 9.903 eV and h-^ = 1.533 eV. AiiJ^ and 
AE'tt are defined as the difference with the bond energy between the structure with a vacancy and 
the perfect structure. A£^o.+7r represents ^E^ + /^.E-,^. 



Total Ideal FN-SN SN-TN SN-TN' Others 







0.9116 


0.9161 


0.9058 


0.9107 








0.1036 


0.1412 


0.1034 


0.1020 




£;^(eV/bond) 




-18.055 


-18.143 


-17.941 


-18.036 




£;^(eV/bond) 




-0.635 


-0.866 


-0.634 


-0.625 




Number of bonds 


124 


124 


12 


12 


24 


76 


^a(cV) 


-2237.918 


-2238.717 


-217.718 


-215.288 


-432.866 


-1372.046 


E^{eY) 


-81.291 


-78.740 


-10.386 


-7.605 


-15.000 


-48.300 




-2319.209 


-2317.457 


-228.105 


-222.893 


-447.866 


-1420.362 


AE^ 


0.799 




-1.068 


1.362 


0.434 


0.071 


AE^ 


-2.551 




-2.766 


0.015 


0.240 


-0.040 


AEcr+TT 


-1.752 




-3.835 


1.377 


0.674 


-0.032 


Ai^.tlTotal) X 100(%) 


-45.6 




61.0 


-77.7 


-24.8 


-4.1 


Ai^.ttTotal) X 100(%) 


145.6 




157.9 


-0.9 


-13.7 


2.3 


Ailgfci) X mm 


100 




218.9 


-78.C) 


-38.5 


1.8 
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